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We study the dynamic behavior of the Ziff-Gulari-Barshad (ZGB) irreversible surface-reaction 
model around its kinetic second-order phase transition, using both epidemic and poisoning-time 
analyses. We find that the critical point is given by p\ — 0.387 368 2 ± 0.000 001 5, which is lower 
than the previous value. We also obtain precise values of the dynamical critical exponents z, 5, and 
rj which provide further numerical evidence that this transition is in the same universality class as 
directed percolation. 
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There has been a great deal of interest surrounding the critical behavior of non-equilibrium kinetic models, including 
directed percolation (DP) the contact process [| and various surface-reaction or catalysis models (for review see 
Q). These models all contain a similar continuous (second-order) phase transition from an adsorbing to a vacuum 
7-H state |,§. 

It has been shown that many of these transitions behave in a universal manner even though the systems abide by 
different local rules and are inherently modeling different physical systems. Grassberger and Janssen || postulated 
that all single-component continuous transitions fall into the robust DP or "Reggeon field theory" class, and many 
numerical simulations have supported this hypothesis (i.e., ||f7|]9|). Grinstein et al. ]5| were the first to hypothesize 
■ that the specific oxygen-poisoning (second-order) transition of the ZGB surface reaction model jll] falls into this class 
| of models, and Jensen et al. ran simulations which support this conclusion [Q. In this communication, we report on 
C^) i new, very extensive simulations which provide further support for this hypothesis, and correct an apparent error in 
the reported value of the location of that transition. Assuming the identification of the ZGB model with the DP class 
to be exact, our results give the most accurate values of the DP dynamical critical exponents to date. 

The ZGB model is a simplified model for the irreversible reaction of CO (A) and O2 (B2) catalytic reaction on a 
Pt surface. The simulation involves the adsorption and reaction of species on a square lattice and proceeds via the 
Langmuir-Hinchelwood mechanism, in which all molecules must adsorb before they can react. The following kinetic 
\ scheme is employed, 

A + *^A* (1) 
B 2 + * -> 2B* (2) 
A* + B* —> AB + 2 * (3) 



where * refers to a lattice site. A Monte-Carlo algorithm is employed where a site is randomly chosen. If the site 
is empty, an A will adsorb with probability p. With probability 1 — p, a B2 is adsorbed and instantly dissociates 
onto that site and a randomly chosen neighboring site if the latter is empty. When a species adsorbs, it checks for 
adjacent neighbors of the opposite species. If one is present, the two react immediately, implying an infinite reaction 
rate as compared to the adsorption rate. There are two transition points in this system. At p 2 , there is a first-order 
(discontinuous) transition to an A-poisoned (saturated) state, and at p\ there is a second-order (continuous) transition 
to a B-poisoned state. Between these points exists a reactive window where the system can reach a steady state and 
react indefinitely. (For a phase diagram, see pd|.) Even within the window, for finite systems, the system is only 
metastable as it can in principle poison by a statistical fluctuation to a non-reactive state. However, here the average 
time to poison t p grows exponentially with lattice size L, the signature of a reactive steady state Jl^j . The value of p 2 
has been accurately determined to be 0.525 60 ± 0.000 01 [[l3| using the constant-coverage ensemble algorithm. This 
algorithm, however, is only applicable to finding the location of the first-order transition. 

Because the second-order transition is a continuous one to a single adsorbing state, it is expected to fall into the 
DP class ||. Indeed, while the ZGB model involves three components (A, B, and vacant sites), at the second-order 
transition, there are rarely A molecules at the surface so it is essentially a two-species model like other members of the 
DP class. The value of its transition point p\ was first empirically observed to be 0.389 ± 0.005 jll|]. A more precise 
value 0.390 65 ± 0.00010 was obtained by Jensen et al. by using an epidemic analysis |?]]. However, while recently 
performing some other investigations ]13[] , we found that this value appears to be somewhat high. Thus, we carried 
out new simulations, using the epidemic procedure as well as a poisoning-time analysis, to re-examine the value of p\ 
and the related dynamic critical exponents. 
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The epidemic method was initially used to study the contact process || and has been successfully applied to 
determine the critical exponents and the critical point for DP [jl[. To run the epidemic analysis, we started with a 
large (1024 x 1024) system completely saturated with B except for a single vacant site in the center. A large system 
is necessary so that the cluster never hits the boundary. The simulation was run at a set value of p and a reactive 
cluster was grown and watched until the system reverted to a non-reactive adsorbate state, or a maximum cutoff time 
was reached. 

The vacant sites (numbering n v ) were kept on a list which was randomly accessed for each adsorption trial, in- 
crementing the time f by l/n v . As each cluster grew, the quantities of interest were recorded in log2 bins of time. 
Since only approximately 3% of all clusters grown reached the last bin, it was necessary to make numerous runs order 
to obtain satisfactory statistics. For the values p = 0.387 36 and 0.38737, 8 x 10 7 clusters (N) were grown up to 
2 13 = 8192 time-steps, requiring a total of 200 days of computational time on a HP 9000 series UNIX platform. In 
the work of j?J in contrast, only 100 000 to 250 000 clusters were grown up to 1000 timesteps. Although we could 
pinpoint p\ to four significant figures in just a few hours, we decided to carry out extensive runs in order to find p\ 
to six significant digits and to determine the dynamical critical exponents precisely. 

We measured the three quantities introduced by Grassberger and de la Torre ||: the survival probability P(t), the 
mean number of vacancies (averaged over N) n(f), and the mean-square radius of gyration of vacant sites (averaged 
over N alive at f), R 2 (t). At the critical point, these are hypothesized to follow the asymptotic power laws, 

P(t)~t~ s (4) 
n(t)~ i" (5) 
R 2 (t)~t z (6) 



These exponents follow the hyperscaling relation 



dz = 2r] + AS (7) 



where d is the spatial dimension. These relationships provide a powerful method to determine p\ by evaluating the 
effects of slightly non-critical values of p, in which case the resulting behavior deviates from ([|) - (g) for large t. An 
example of this is shown in Fig. 1 where n(f) is plotted for p = 0.390 65, 0.387368 2, and 0.38407. The upper curve 
is for the value of p\ reported in 0, while the central curve is for the value we find below. 
In order to find the exponents accurately, we consider the local slopes which are defined as 

-5(f) =ln[P(f)/P(f/2)]/m2 (8) 

and similarly for 77(f) and z(t). (Here we used a factor of 2 rather than 5 or 8 of previous work 0,f7j, which we could 
do because of our higher statistics.) These are all graphed in Fig. 2 for p = 0.387 36 and 0.387 37. The local slopes 
can be expanded as Q| 

S(t)=S+j + ^r. (9) 

If the non-analytic corrections were negligible, then it would be easy to extrapolate the critical exponents as a function 
of 1/f as discussed by Grassberger. However, these corrections are rather large and therefore hinder a direct linear 
extrapolation. In order to overcome this problem, we grew over 5 • 10 8 clusters to 2 9 time-steps so that we could 
better follow the non-analytical trajectory of each curve. Extrapolating these results to f — > 00, we find 

6 = 0.4505 ± 0.001, 77 = 0.2295 ± 0.001, z = 1.1325 ± 0.001. (10) 

consistent with 5 = 0.452 ± 0.008, 7/ = 0.224 ± 0.10 and z = 1.133 ± 0.002 found in @. 

For comparison, the updated values recently found by Grassberger and Zhang || for DP are 

5 = 0.451 ± 0.003, 77 = 0.229 ± 0.003, z = 1.133 ± 0.002. (11) 

The precise agreement between ([!(]) and ( fil| ) leaves little doubt that the ZGB model is included in the DP class as 
predicted by Grinstein et al. [|| . 

For p away from pi, 71(f) follows the scaling behavior || 

77(f) -f-VKp-Pl)^"] (12) 

and similarly for P(f) and R 2 (t). It follows from this equation that 
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d\np 



oc t 1 '^ (13) 

p=pi 



where vu is the time-domain correlation length exponent. In Fig. 3, we plot the quantity on the left-hand side of the 
equation above, calculated by taking the difference of n(t) for p = 0.387 36 and 0.38737, vs. Int. This plot shows that 
Grassberger's value for DP v\\ = 1.295 ± 0.006 Q is completely consistent with our data. 
To determine the precise value of p\ , we expand the scaling function (f> as 

Pit) «^[a + 6(p- Pl )i 1/l/ ii ...] (14) 

This equation implies that a plot of P(t)t s vs. t}-l v W for values p close to p\ should yield straight lines and that a direct 
linear interpolation of the data from different values of p can be used to find p\ (which corresponds to a horizontal 
line in this plot) (Fig. 4). There is an initial curvature which is to be expected for small clusters due to finite-cluster 
effects. To minimize this effect, we have added a constant c to the time that effectively allows for a higher-order 
analytic correction term: 

it + c)- s *t- s (l- S -f) (15) 

where c « 1.7 was found to give the best results. The resulting plot of our data is shown in Fig. 4. The statistical 
fluctuations in each bin are given by 



V AWi (16) 

which implies that the largest bins which have the most accurate data also have the greatest error (least precision). 
Interpolating the two data curves in Fig. 4, we deduce that p\ is given by 

Pi = 0.387 368 2 ±0.000 001 5. (17) 

This result is nearly two orders of magnitude more precise than the result of @, 0.390 65 ± 0.000 10, and more than 
30 combined error bars lower. We believe that some error must have occurred in the simulations or analysis of 0). 

To confirm our value for pi, we also ran a poisoning-time analysis |l3| of the system at its critical point. Similar 
methods have been applied to other problems including the quantification of finite lattice effects JT^,^,Q. To do this, 
we essentially run the opposite dynamic algorithm performed by the epidemic analysis. We start with a small lattice 
in a fully reactive state (all vacant sites) and set the value of p at our determined p\. Periodic boundary conditions 
are applied and the system is allowed to run until the adsorbate B saturates or poisons the system, causing a global 
non-reactive state. When the value of p is at pi, it is expected that the dependence of t p on L will be power-law, and 
when p ^ pi, the dependence will be exponential fl^ , |l4}| . We ran this simulation for square lattices of powers of 2 in 
size from 8 x 8 to 64 x 64 for roughly 10 5 runs each. Here, a time-step is defined as L 2 adsorption trials. Fig. 5 shows 
the results of our analysis and it was found that at p\ , the relationship is indeed 



t P ~ L w (18) 

with w — 1.77±0.02. In we observed that w — 2/z — v\\/v^, indicating that the time to expand a reactive state 
scales as the time to contract. The z implied by this result is consistent with the value determined above. While this 
method is evidently less efficient than the epidemic analysis, it provides a useful confirmation our results for p\ and 
z. 

In conclusion, we have provided improved numerical evidence that the ZGB oxygen-poisoning transition falls into 
the larger DP class of nonequilibrium models. Accepting that that hypothesis is true, which seems certain, our 
exponents represent the most accurate values of the DP dynamic critical exponents to date (by a factor of about two) . 
We also independently confirm that the value of for the ZGB model falls into the DP class and use it to find a 
highly accurate and corrected value of p\ . 

This material is based upon work supported by the US National Science Foundation under Grant No. DMR-9520700. 
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FIG. 1. The behavior of the number of vacant sites n as plotted against time t, for p — 0.390 65, 0.387 368 2, and 0.384 07 
(top to bottom). 

The upper value is pi from M, and the center is for the value found here. 



FIG. 2. The three critical exponents derived from our epidemic analysis: S (a), r/ (b), and z (c). These values show 
super-critical (o: p = 0.387 37) and sub-critical (A: p — 0.387 36) behavior. Each of these lines represents the average of 
3.5 x 10 7 runs. The actual value of pi falls between these lines and can be determined by linear interpolation as in Fig. 4. 

FIG. 3. Plot to determine u« from (13) by using the values of p = 0.387 36 and 0.387 37. The line is for vu — 1.295 as 
determined by Grassberger for DP. It can be seen that the ZGB data is consistent with this value. 

FIG. 4. Plot allowing a linear interpolation for p\ as expressed in equation (14). The lines for 0.387 36 (A) and 0.387 37 
(o) represent the sub- and super-critical behavior respectively. The bold line represents the interpolation for the value, 
Pi — 0.387 368 2. The error bars were calculated as in (16). Here, t is offset by an additive constant 1.7 to improve small-time 
behavior. 

FIG. 5. Results of our poisoning-time analysis for the same values of p displayed in Fig. 1 (A: 0.390 65, o: 0.387 368 2, o: 
0.384 07). This plot demonstrates that the expected power-law behavior obtains when our value of pi is used. 
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